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I.  INTRODUCTION  AND  SUMMARY 


We  have  assembled  a  set  of  three  computer  codes  and  interfaced 
them  to  each  other  in  order  to  retrieve  approximate  surface  wave 
path  corrections.  They  are  (1)  TELVEL,  an  interactive  time  series 
analysis  program  which  extracts  dispersion  curves  and  spectral 
amplitudes  from  observed  wave  trains,  (2)  SWIP,  which  inverts 
dispersion  curves  and  attenuation  function  in  terms  of  an  average 
earth  structure  along  the  path,  and  (3)  SURWAV  which  generates 
synthetic  surface  wave  trains  and  spectra  for  comparison  with 
observations. 

The  analytical  procedure  for  constructing  path  corrections 
comprises  four  major  steps: 

1.  Group  and  phase  velocity  dispersion  curves  are  extracted 

from  each  seismogram.  A  narrow  band  filtering  technique 
yields  an  initial  estimate  of  group  velocity  as  a 
function  of  frequency.  This  estimate  is  then  refined 
through  an  iterative  technique  based  on  phase  matched 
filters  (Herrin  and  Goforth,  1977),  which  permits 

retrieval  of  group  and  phase  velocity  dispersion  curves 
after  removal  of  the  most  severe  contamination  due  to 
multipathing. 

2.  The  dispersion  curves  are  then  inverted  for  earth 
structure  along  the  path,  as  described  by  Bache,  Rodi  and 
Harkrider  (1978).  Since  the  data  are  most  sensitive  to 
the  shear  velocity  structure,  density  and  compressional 
velocity  are  usually  constrained  by  an  empirical  relation 
such  as  Birch's  law.  The  tradeoff  curve  between  model 
norm  and  r.m.s.  data  misfit  can  be  outlined  by  performing 
a  single  iteration  for  each  value  of  the  tradeoff 
parameter  e  and  successive  iterations  are  performed  for 
the  optimal  value  of  e  until  convergence  is  attained. 
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3.  We  use  this  inversion  model  to  compute  theoretical 
unnormalized  spectral  amplitudes.  In  the  absence  of 
strong  lateral  variations,  these  theoretical  amplitudes 
should  differ  from  the  observed  ones  only  by  a  constant 
proportionality  factor  times  an  effective  attenuation 
factor  of  the  form  exp[-Y(«)r]. 

4.  Comparison  of  observed  and  synthetic  spectra  then  permits 
us  to  estimate  simultaneously  the  source  strength  from 
the  spectral  ratio  at  long  periods,  and  the  effective 
attenuation  factor  as  a  function  of  frequency,  which  can 
be  interpreted  in  terms  of  an  average  Q  structure  for  the 
path.  The  overall  effectiveness  of  the  procedure  can 
then  be  assessed  by  comparing  observed  and  synthetic 
seismograms. 

The  results  from  steps  1  to  3  are  relatively  robust  with 
respect  to  minor  errors  and  noise  in  the  data  and  are  thus  reliable 
in  most  cases.  On  the  other  hand,  step  4  is  much  more  sensitive  to 
lateral  structural  variations  and  unresolved  multipathing  which  are 
not  accounted  for  in  the  synthetic  calculations,  so  that  Q  estimates 
are  much  more  questionable.  The  simplest  approach  for  step  4  uses 
the  fact  that  y(<d)  should  vanish  in  the  long  period  limit.  We  have 
verified  that  this  constraint  provides  acceptable  estimates  of 
moment  and  attenuation  function  for  low-noise  records. 

We  have  tested  this  sequence  of  programs  on  a  variety  of 
synthetic  seismograms,  and  on  two  real  seismograms. 

Synthetic  records  involving  simulated  multipathing  were 
constructed  using  a  monopole  (explosion)  source  and  an  attenuating 
earth  model  valid  for  the  eastern  U.S.  Multipathing  was  simulated 
by  superposition  of  synthetic  wave  trains  calculated  for  several 
epicentral  distances.  We  applied  the  procedure  outlined  above  to 
these  synthetic  records  with  the  following  results. 
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1.  Multipath  arrivals  were  successfully  identified  and 
separated  by  the  algorithm  whenever  the  corresponding 
time  delay  was  about  30  seconds  or  more.  More 
importantly,  even  when  adequate  time  resolution  could  not 
be  achieved,  the  ensuing  bias  in  dispersion  estimates  was 
found  to  be  small.  For  extreme  multipathing  (e.g.,  many 
narrowly  spaced  arrivals)  the  procedure  failed  to  yield  a 
recognizable  dispersion  curve,  but  did  not  result  in  a 
stable,  but  misleading  dispersion  cur»e. 

2.  The  inversion  of  group  and  phase  velocity  dispersion 
recreated  the  initial  "true"  structure  quite  accurately. 
Numerical  experiments  show  that  the  selection  of  a 
position  along  the  tradeoff  curve  close  to  the  optimum  is 
fairly  important,  and  that  the  step  involving  approximate 
construction  of  the  tradeoff  curve  should  not  be  bypassed. 

3.  Spectral  amplitudes  of  synthetics  calculated  from  the 

(unattenuating)  inversion  model  match  closely  the 
corresponding  amplitudes  for  the  "true"  structure  when 
attenuation  is  ignored  in  the  latter.  Consequently, 
comparison  of  synthetic  spectra  for  the  inverted 
structure  with  "observed"  spectra  yielded  an  estimate  of 
explosion  moment  accurate  to  -  2  percent,  as  well 

as  reasonably  accurate  estimates  of  y(id)  at  all 
frequencies  (Table  I). 

Thus,  in  the  absence  of  complications  other  than  those 
included  in  this  test,  the  method  appears  to  yield  acceptable  path 
correction  estimates.  As  a  considerably  more  severe  test  of  the 
procedure,  we  also  analyzed  two  seismograms  from  the  Hoggar 
explosion  SAPHIRE,  recorded  at  AAE  and  SHI. 

We  were  able  to  retrieve  stable  dispersion  curve  estimates  in 
the  range  10-50  seconds  from  both  records,  and  to  invert  those 
curves  for  average  path  structures.  However,  comparison  of  observed 
and  synthetic  spectral  amplitudes  reveals  significant  contamination 
from  effects  not  included  in  the  forward  modeling.  Specifically, 
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the  AAE  seismogram  seems  to  be  deficient  in  long  period  energy,  in 
contrast  to  the  SHI  seismogram  which  has  excess  long  period  energy. 
The  SHI  record  exhibits  long  period  drift  whicn  may  reflect  some 
instrumental  problem.  In  the  case  of  AAE,  we  speculate  that  the 
proximity  of  this  station  to  the  African  Rift  Zone  may  be  a  source 
of  bias  which  is  not  handled  well  by  the  procedure. 

Although  we  can  formally  estimate  y{m)  and  u'J'oo  from  such 
data,  it  is  clear  that  these  estimates  cannot  be  very  reliable,  and 
that  a  more  constrained  approach  is  called  for  whenever  such 
difficulties  arise.  Interestingly,  such  formal  estimates  lead  to  a 
fair  match  between  observed  and  synthetic  wave  trains.  We  are 
considering  two  alternatives.  The  first  one  is  to  restrict  the 
range  of  allowable  Q  structures  based  on  independent  geophysical 
information,  and  to  achieve  a  best  fit  to  the  observations  within 
this  allowable  range  (e.g.,  Cheng  and  Mitchell,  1981).  Another 
alternative  consists  of  a  simultaneous  interpretation  of 
observations  from  multiple  event-station  pairs  in  terms  of  a 
regionalized  Q  model  of  the  crust  and  upper  mantle.  In  either  case, 
additional  information  from  independent  evidence  is  called  for. 

Current  efforts  are  aimed  at  simplifying  input,  documenting 
and  streamlining  the  interfacing  of  these  three  codes  so  as  to 
permit  their  transfer  to  the  VSC  seismic  system  and  begin  further 
testing  using  high  quality  digital  data. 
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II.  THEORETICAL  SURFACE  WAVE  AMPLITUDES 


Harkrider  (1964,  1970)  formally  solved  the  problem  of 
computing  synthetic  seismograms  from  a  buried  source  in  a 
plane-layered  earth.  Comparison  of  observed  seismograms  with 
synthetic  seismograms  allows  the  formulation  of  the  inverse  problem 
of  obtaining  earth  structure,  dispersion  and  source 
characteristics.  The  vertical  component  of  the  Rayleigh  wave 
displacement  (up  »  positive)  from  an  explosion  at  depth  h  is  given 
by  the  expression  (Bache,  Rodi  and  Harkrider,  1978): 


W(r,«) 


-4irUj  ^(iu) 


Ks(h)  arU) 
c 


where 

w 

r 

0) 

us 

'■P 

Ks(h) 

Ar(u) 

C 

«<2) 

y(ui) 

Re 

A 


Vertical -component  of  displacement  for  funda¬ 
mental  mode  Rayleigh  wave. 

Event-station  separation  (km). 

Angular  frequency. 

Source  layer  shear  modulus. 

Reduced  velocity  potential. 

Source-depth-dependent  modal  excitation  factor. 
Path  structure  response. 

Phase  velocity. 

Hankel  function  of  the  second  kind  order  zero. 
Anelastic  attenuation  factor. 

Earth  radius. 

Event  station  angular  separation. 


All  factors  in  this  expression  are  real  and  positive  except 
for  the  Hankel  function  and  P{u) .  All  phase  information  is 
therefore  contained  in  these  two  factors.  With  >>  1,  the 
asymptotic  approximation  of  the  Hankel  function  is 
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Furthermore,  for  small  values  of  w,  the  pha^e  of  t j)(u»)  is  small;  so 
almost  all  of  the  phase  information  is  contained  in  the  exponential 
function.  Given  the  distance  r,  the  phase  velocity  c(<u)  and  group 
velocity  u(«)  are  recovered  from  the  phase  spectrum  using  the  method 
of  phase-matched  filters  (Section  III). 

The  phase  and  group  velocity  dispersion  curves  are  inverted 
for  earth  structure  following  the  method  of  Bache,  Rodi  and 
Harkrider  (1978).  This  inversion  depends  only  on  the  phase  function 
of  the  seismogram  (Figure  1).  The  earth  structure  may  then  be  used 
to  compute  the  unknown  quantities  AR(w)  and  Ks(h)  in  Equation 
(1).  Of  course,  the  source  depth  h  is  usually  not  known,  but  since 
K$(h)  is  a  slowly  varying  function  for  small  h,  an  approximate 
depth  estimate  is  usually  adequate. 

In  the  long  period  limit  <!»(«■>)  38  a  positive  constant. 
Evaluation  of  Equation  (1)  at  an  observation  point  with  y(w)  *  0  and 
'I'oo*  1  allows  us  to  form  the  ratio 


usl*o°c-Y(M)r 

uinv 


(3) 


where  u-„.,  is  the  value  of  the  shear  modulus  u  for  the  source 
inv 

layer  obtained  from  the  inversion.  If  the  records  are  of  high 
quality,  then  Equation  (3)  directly  determines  us>^oo  and  y(w)  by 
using  the  fact  that  y(«)  *  0« 

Finally,  the  attenuation  coefficients  y(<*>)  can  be  inverted  to 

determine  Q  as  a  function  of  depth,  at  least  in  principle.  In  terms 

of  the  quality  factors  Q  and  Q.  of  individual  layers  in  the 

a  6 

model,  the  attenuation  coefficients  may  be  written  (Anderson,  et 
al.,  1965) 
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Figure  1.  Major  factors  influencing  surface  wave  spectra.  The  amplitude  and  phase  factors  are 

divided  Into  source-  and  path-dependent  features.  A  full  path  correction  must  determine 
the  three  features  on  the  right  of  the  diagram. 


where  the  index  l  represents  the  layer  number.  The  derivatives  in 
Equation  (4)  are  determined  in  the  inverse  calculation  for  earth 
structure.  The  phase  velocity  is  largely  insensitive  to  variations 
in  a;  so  the  second  sum  is  small  compared  to  the  first  one.  In 
addition,  if  all  losses  occur  in  shear,  with  no  bulk  attenuation, 

Q 

then  Q  5s  tQ„»  which  can  be  used  as  a  constraint.  Failure 
to  obtain  a  physically  reasonable  model  of  attenuation  structure  is 
usually  diagnostic  of  a  complication  in  wave  propagation  (e.g., 
multipathing,  lateral  inhomogeneities)  or  in  the  recording  (e.g., 
noise,  faulty  instrumentation),  which  has  not  been  correctly 
identified  in  the  analysis. 

For  seismograms  which  are  too  contaminated  by  noise  to  use 
Equation  (3)  for  estimating  y(«)»  we  may  still  take  attenuation  into 
account,  but  the  range  of  allowable  Q  models  must  be  severely 
restricted,  at  the  cost  of  achieving  a  less  satisfactory  fit  to  the 
observed  seismogram.  This  procedure  was  adopted  by  Cheng  and 
Mitchell  (1981),  who  consider  a  one  parameter  family  of  Q  models  in 
their  simultaneous  analysis  of  several  events  with  known  moments. 

The  combined  forward  and  inverse  modeling  used  to  obtain  the 
path  corrections  is  summarized  in  Figure  2  in  the  form  of  a  flow 
diagram. 
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FORWARD  MODELING 

OBSERVED  SEISMOGRAM  CONSISTENCY  CHECK  — SYNTHETIC  SEISMOGRAM 


Figure  2.  Flow  diagram  for  determining  path  corrections.  Dispersion  curves  and  amplitudes  are  found 
by  phase-matched  filtering.  Synthetic  amplitudes  from  inverse  structure  are  compared  with 
observed  amplitudes  to  estimate  source  amplitude  and  attenuation.  Attenuation  coefficients 
are  constrained  by  and/or  inverted  for  a  Q  model.  The  output  is  a  geophysical  model  (velo¬ 
city  and  Q  structure)  and  the  observed  path  effects  (dispersion  and  attenuation). 


III.  PHASE-MATCHED  FILTERING 


The  technique  of  phase-matched  filtering  for  obtaining 
surface-wave  dispersion  and  spectral  amplitudes  was  first  proposed 
by  Herrin  and  Goforth  (1977),  although  Dziewonski,  et  al^  (1972) 
used  essentially  the  same  technique  without  iteration.  Earlier 
methods  included  the  graphical  "peak  and  trough"  method,  narrow-band 
filtering  (for  example,  Archambeau,  Flinn  and  Lambert,  1966; 
Dziewonski,  Block  and  Landisman,  1969),  and  moving-window  analysis 
(Landisman,  et  al.  1968;  Landisman,  Dziewonski  and  Sato,  1969).  Of 
these,  the  peak  and  trough  method  is  simple,  but  rapidly  becomes 
unreliable  in  the  face  of  complications  in  the  seismograms. 
Dziewonski,  et  al .  (1972)  demonstrate  that  the  other  two  methods 
lead  to  systematic  errors  in  the  dispersion  estimate  whenever:  (a) 
the  amplitude  spectrum  A(w)  varies  with  frequency,  and  (b)  the 
phase  velocity  c(«)  differs  from  the  group  velocity  u(w).  Rarely, 
if  ever,  are  these  two  conditions  not  satisfied. 

The  Herrin  and  Goforth  method  is  based  upon  classical 
matched-filter  principles.  The  basic  idea  is  to  remove  as  much 
dispersion  from  the  seismogram  as  possible  so  that  the  resulting 
impulse-like  signal  becomes  easy  to  process  with  conventional 
techniques. 

Suppose  a  fundamental -mode  surface  wave  s(t)  has  a  Fourier 
transform 


S («)  -  |S(u)|  elk(w)r  . 

The  phase  and  group  velocities  are 


cU) 


i  \  du 

U(u)  "  TOTF 


(5) 

(6) 
(7) 
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We  seek  an  estimate  of  k(u»).  Suppose  a  known  dispersion  curve  is 
specified  by  k‘ («): 

S 1  (a.)  -  |S'(u)|eikl(w)r.  (8) 

The  cross-spectrum  of  S'(w)  and  S(<u)  is  given  by 

H(ci»)  -  S(u>)  S'*(«0  -  |S («) j  |S'(«)|e^k{“}  '  k'(u);ir  (9) 

where  *  denotes  complex  cisjugation.  Ideally,  if  k'(«)  »  k(w),  the 
cross-correlation  is  phaseless  (an  even  function).  In  addition,  if 
|S'(«)|  *  l/|S(«)l,  then  |h(«)|  ■  1  and  h(t)  is  an  impulse.  S(u>)  is 
then  the  transfer  function  of  a  deconvolution  or  inverse  filter. 
Alternatively,  if  S'(w)  *  S(w),  H(w)  *  |s(u)|^  and  h(t)  is  the 
autocorrelation  function  of  s(t).  The  salient  feature  to  note  is 
that  H(u)  has  much  less  phase  than  S(w),  so  that  h(t)  is  compressed 
in  time  with  respect  to  s(t). 

The  residual  phase  in  H(«)  is  small  if  k'(<i>)  s=k(u),  and 
consequently,  the  phase  spectrum  of  H(w)  is  much  easier  to  unwrap 
compared  to  the  phase  spectrum  of  S(u»).  If  0h(u>)  is  the  phase 
spectrum  of  H(u)),  then  k(w)  *  k‘(w)  +  >)/r. 

Because  k'(<u)  is  chosen  to  approximate  k(n>),  S'(ui)  is  called  a 
phase-matched  filter  (PMF).  The  amplitude  spectrum  of  S'(u>)  can  be 
arbitrarily  chosen.  The  deconvolution  filter  |S'(w)|  *  |S( ui~)  1 
has  good  time  resolution,  but  is  extremely  sensitive  to  extraneous 
noise.  When  |S ' (u)  |  *  |S(w)|,  then  S'*U)  is  the  classic  matched 
filter  which  optimizes  the  signal-to-noise  ratio  in  a  least-squares 
sense.  The  resulting  time  resolution,  however,  is  poor.  A  useful 
compromise  between  these  extremes  which  the  code  utilizes  is 
|S'(«)|  *  1  so  that  |H(o)>  |  -  |S(»)  1.  Since  the  cross-spectrum 
preserves  the  amplitude  spectrum  of  the  observed  seismogram,  the  PMF 
output  spectrum  yields  noise-suppressed  spectral  amplitude  estimates. 

The  greatest  advantage  of  the  phase-matched  filtering  approach 
stems  from  the  time  compression;  it  affords  better  noise  rejection 
in  spectral  amplitude  estimation  and  multipath  separation.  Both 
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these  topics  are  well  covered  by  Herrin  and  Goforth  (1977),  and  will 
be  only  briefly  mentioned  here. 


Since  the  filter  is  matched  to  the  seismogram  and  not  to  the 
noise,  the  signal  energy  is  compressed  to  times  near  t  ■  0  while  the 
noise  is  not.  Therefore,  a  spectral  estimate  of  h(t),  evaluated 
from  a  time  window  near  t  -  0,  is  less  contaminated  by  noise.. 

A  multipathed,  or  a  higher-mode  arrival  with  group  arrival 
time  tM  will  be  time-shifted  to  t  ■  tM  -  tp  where  tp  is  the 
primary  group  time,  and  not  to  t  »  0.  In  addition,  if  the 
dispersion  curve  of  the  multipath  does  not  match  that  of  the 
primary,  then  the  multipath  component  of  the  output  signal  will  not 
be  symmetrical.  Herrin  and  Goforth  give  illustrative  examples  of 
these  two  effects. 

Of  critical  concern  is  the  resolving  power  of  the  PMF  method; 
namely,  how  much  time  separation  between  primary  and  multipath 
arrivals  is  necessary  to  yield  unperturbed  dispersion  curves  and 
spectral  amplitudes  for  the  primary.  The  time  width  of  the  output 
h(t)  depends  on  the  bandwidth  of  |h(w)|,  which  in  turn  depends  upon 
the  seismogram  bandwidth.  A  larger  bandwidth  gives  a  narrower  time 
domain  signal,  according  to  the  uncertainty  relation  at  a<o  ==  1.  Our 
experience  with  synthetic  data  shows  this  to  be  the  case. 

The  block  diagram  for  the  code  TELVEL  is  shown  in  Figure 
3.  The  algorithm  is  similar  to  that  of  Herrin  and  Goforth,  except 
that  the  initial  group  velocity  estimate  is  obtained  by  narrow-band 
filtering  instead  of  being  user-provided. 
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Figure  3.  Block  diagram,  dispersion  curve  extraction.  The  program  TELVEL  uses  phase-matched  filtering 
as  proposed  by  Herrin  and  Goforth  (1977). 


IV.  TESTS  OF  PHASE-MATCHED  FILTERING  USING 
SYNTHETIC  SEISMOGRAMS 


Our  tests  with  synthetic  data  started  with  the  model  and 
associated  dispersion  curves  shown  in  Figure  4.  Seismograms  were 
generated  for  various  epicentral  distances  and  additively  superposed 
to  simulate  multipathing  (Figure  5).  The  source  was  an  explosion  at 
800  meter  depth  with  ^00»  104  m^. 

As  expected,  TELVEL  performed  very  well  on  synthetic  SI,  which 
involves  no  multipathing,  and  on  synthetic  S2,  which  has 
well-separated  groups.  Typical  phase  velocity  errors  were  within 
^  0.002  km/sec,  and  typical  group  velocity  errors  within 
+  0.004  km/ sec  in  the  5-200  sec  band.  Greater  group  velocity  errors 
are  to  be  expected  because  of  the  derivative  operation  (Equation 
7).  Figure  6  shows  the  clean  separation  of  arrivals  achieved  for 
synthetic  S2. 

Test  case  S3  is  at  the  limit  of  resolution  of  the  phase 
matched  filter  method.  The  first  two  arrivals  could  not  be  clearly 
separated.  Nevertheless,  S3  yielded  good  phase  velocities  (errors 
<  +  0.005  km/sec)  in  a  limited  period  range  of  5.3  to  40.0  seconds. 
Group  velocity  errors  were  about  +  0.010  km/sec  from  5.7  to  25.0 
seconds. 

Interference  of  the  closely  spaced  arrivals  in  S4  made  it 
impossible  to  obtain  a  dispersion  curve  at  all  for  this  seismogram. 
Individual  arrivals  could  not  be  resolved  using  either  phase-matched 
filtering,  or  narrow-band  filtering.  Even  in  this  extreme  case, 
however,  the  complexity  of  the  seismogram  did  not  lead  to  an 
incorrect  dispersion  curve;  it  was  simply  impossible  to  obtain  any 
dispersion  curve. 

Synthetic  S5  was  the  only  one  in  this  series  where  higher 
modes  were  included.  The  analysis  of  this  seismogram  presented  no 
difficulty  —  probably  because  the  higher  modes  were  of  low 
amplitude  and  distinct  from  the  primary  arrival  in  both  frequency 
and  time.  Errors  were  comparable  to  those  obtained  with  SI  and  S2. 
These  tests  on  synthetic  seismograms  show  that  the  PMF  method  is  an 
accurate  and  stable  way  of  finding  dispersion  curves. 
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Figure  4.  Eastern  United  States  structure  and  Rayleigh  fundamental -mode  dispersion  waves 
structure  Is  model  SI  of  Bache,  Swanger  and  Skholler  (1980). 
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PHF  output  for  S2.  The  two  sharp  pulses  represent  the  two  separate  arrivals.  The  second 
pulse  Is  asymmetrical  because  the  PMF  Is  matched  only  to  the  dispersion  of  the  primary 
arrival . 


It  should  be  noted  in  the  foregoing  discussion  that 
single-station  phase  measurements  are  ambiguous  with  respect  to 
additive  multiples  of  +  2irn,  since 

ls(„)|eikr-  |s(.)|eiCk-lr-:|r  .  (10) 

The  correct  value  of  n  can  be  found  if: 

1.  The  phase  velocity  at  one  frequency  is  known; 

2.  The  data  extends  to  low  enough  frequency  so  that  the  fact 
that  ^  k(u)r  »  0  can  be  exploited; 

3.  Multiple  station  data  of  proper  spacing  and  bandwidth  are 
available. 

An  improper  value  of  n  affects  the  phase  velocity  estimates 
(particularly  at  low  frequencies  since  c  *  ^  +rgnff-)  but 
does  not  effect  group  velocity  estimates  (Equation  7). 
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V  STRUCTURE  INVERSION  -  TEST  ON  SYNTHETIC  DATA 


The  observed  dispersion  curves  can  be  systematically  inverted 
to  infer  the  average  earth  structure  along  the  path  of  travel. 
Although  the  relation  between  dispersion  and  structure  is  nonlinear, 
the  problem  can  be  linearized  and  a  best-fit  solution  to  the  data 
approached  by  successive  iterations.  This  method  was  used 
successfully  by  Bache,  et  al .  (1978)  to  infer  crustal  structure  in 
the  southwestern  United  States  .  The  method  may  also  be  adapted  to 
infer  Qfl(z). 

The  procedure  assumes  that  the  entire  travel  path  from  source 
to  receiver  can  be  adequately  represented  by  a  plane  layered  earth 
structure.  It  yields  estimates  «(z),  s(z),  and  p(z)  corresponding 
to  average  values  of  these  quantities  over  the  path.  In  addition, 
the  corresponding  resolving  kernels  are  computed  and  characterized 
by  a  simple  measure  of  spread.  In  practice,  c(u>)  and  u(n>)  do  not 
permit  estimation  of  a(z)  and  p(z)  independently  of  s(z).  There 
are,  however,  good  geophysical  arguments  for  constraining  a  and  p  to 
be  linear  functions  of  a: 

o(z)  =  Cla  (z)  6(z)  +  C2a  (11a) 

p(z)  =  Clp (z)  a(z)  +  C2(j  .  (lib) 

In  other  words,  we  assume  that  a  Birch-type  law  holds  in  each  layer, 
and  fix  Poisson's  ratio  at  each  depth.  Dobrin  (1976,  page  53)  shows 
empirical  evidence  supporting  these  assumptions. 

The  basic  requirement  imposed  on  a  model  is  that  it  should  fit 
the  data.  However,  measurement  errors,  as  well  as  the  finiteness  of 
the  data  set,  lead  to  well-known  difficulties  and  often  to  unstable 
or  geophysically  implausible  solutions.  In  order  to  ensure 
stability  as  well  as  smoothness  of  the  solution,  the  quantity 

2  2 
(r.m.s.  misfit)  +  9  (model  norm) 
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is,  in  fact,  minimized.  The  model  norm  is  defined  by 


so  that  a  smoother  model  possesses  a  smaller  norm,  e  parameterizes 
the  inevitable  tradeof'  between  these  two  terms,  which  is 

illustrated  on  Figure  7  in  the  form  of  a  "tradeoff  curve.” 

Inversion  theory  as  applied  in  this  report  is  discussed  by 
Wiggins  (1972),  Jackson  (1972),  and  Jordan  (1972). 

The  optimal  point  on  the  tradeoff  curve  is  defined  as  the 
point  for  which  neither  r.m.s.  misfit  nor  model  norm  are  too  large. 
Typically,  one  chooses  a  value  of  e  corresponding  to  the  point  on 
the  curve  closest  to  the  origin,  (e  ~  100  on  Figure  7). 

Other  features  of  our  inversion  program  include: 

1.  The  capability  to  impose  local  constraints  such  as  fixing 
the  velocities  at  discrete  depth  points. 

2.  The  possibility  of  introducing  discontinuities  at 

selected  depths  (e.g.,Moho,  or  sedimentary  layer  if  high 
frequency  data  are  available). 

The  next  step  in  our  synthetic  test  was  to  invert  the 

dispersion  data  obtained  from  SI  (Figure  8).  This  test  was  a 
"blind"  one  in  that  the  operator  did  not  know  the  "real"  structure. 
The  inversion  model  shown  in  Figure  8  incorporated  no 

discontinuities,  and  no  fixed  data  points  were  entered.  The  near 
surface  layers  are  not  resolved  due  to  the  lack  of  short  period 
(T  <  5  second)  dispersion  data. 

The  steep  rise  in  s(z)  in  the  20-45  km  region  is  evidence  of  a 
discontinuity  that  is  smoothed  by  the  inversion.  The  next  inversion 
run  incorporated  a  discontinuity  at  the  midpoint  of  this  region, 
which,  in  fact,  coincided  with  the  "real"  crust/mantle  interface 
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RMS  Misfit  m/sec 


Model  Norm  ||m|| 


Figure  7. 


Tradeoff  curve.  The  RMS  fit  of  the  model  to  the  observed  data 
is  plotted  against  the  model  norm  for  various ivalcJi®Luf 
tradeoff  parameter  0.  For  high  0,  the  model  is  smooth  Ue/3z 
nparl v  constant)  but  the  data  fit  is  poor. 
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Figure  8.  Inversion  model  II,  no  crust/mantle  discontinuity.  The  actual  model  Is  shown  as  a  solid 
line;  the  Inverted  model  as  a  dotted  line.  The  steep  gradient  between  30  and  40  km  Is 
symptomatic  of  a  discontinuity  mollified  by  the  smoothness  criterion. 


(Figure  9).  The  accuracy  of  this  inversion  is  remarkable.  When 
compared  to  the  real  model,  the  inversion  used  a  Birch's  Law  that 
was  significantly  different  from  that  used  to  generate  the  synthetic 
seismograms. 


♦ 


> 


Inversion 
(After  Dobrin) 

Crust  p  »  0.228  +  1.9 
Mantle  p  »  0.65s  +  0.4 


Eastern  United  States 
(Bache  and  Swanger,  1980) 

p  -  0.82s  ♦  0.12 
p  ■  0.65e  +  0.4 


To  fully  illustrate  the  capabilities  of  the  procedure,  we  ran 
a  third  inversion  with  the  "true1'  relations.  In  all  cases, 
Poisson's  ratio  was  fixed  at  0.27  (a  *  1.80s).  This  last 
calculation  gave  a  very  accurate  model,  except  for  the  near-surface 
sedimentary  layers  (Figure  10). 

Although  different  in  structure,  the  three  models  presented 
all  fit  the  data  well.  R.M.S.  misfits  to  the  observed  dispersion 
curves  were  0.006  km/sec  in  all  three  cases  (Figure  11);  so,  on  the 
basis  of  the  data  alone,  the  three  models  are  equally  good. 
However,  without  knowledge  of  the  "true"  structure,  we  would  have 
chosen  the  second  inversion  (Figure  9)  because  it  contains  a 
crust/mantle  structure,  and  it  satisfies  empirical  density-velocity 
relations  discussed  in  the  literature  (e.g.,  Dobrin). 

The  dispersion  data  included  observations  every  0.005  Hz  from 
0.005  Hz  to  0.20  Hz  for  both  c(«)  and  u(w).  The  variance  estimates 
and  resolution  of  the  inversion  model  are  summarized  below  in  terms 
of  standard  deviations  of  the  estimated  parameters  and  spreads  of 
the  resolution  kernels. 


23 


* 


SYSTEMS.  SCIENCE  AND  SOFTWARE 


•g  s>  o>  s’c  o'c  s'z  Q'Z  Q'i  o-i  s-o  cro 


VI 39 


24 


SYSTEMS.  SCIENCE  AND  SOFTWARE 


Figure  9.  Inversion  model  #2,  discontinuity  at  34  km.  On  the  basis  of  Model  #1  (Figure  8) 
discontinuity  was  Inserted  in  the  inversion  model  at  34  km. 


0.0  20.0  40.0  60.0  80.0  100.0  120.0  140.0  160.0  180.0  200.0 

DEPTH 

Figure  10.  Inversion  model  #3,  In  which  the  known  p(0)  relation  Is  used.  This  Inversion  was  run  to 
show  the  capabilities  of  the  Inversion  procedure.  Normally,  the  p(0)  relation  In  a  glvei 
region  Is  not  known,  and  global  average  values  must  be  used  as  In  Model  #2. 


MODEL  PARAMETERS 
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Figure  11.  Dispersion  curves  from  Inversion  model  #2.  These  curves  should  be  compared  to  those  In 
Figure  4.  The  RMS  misfit  Is  0.006  km/sec. 


OEPTH 

(km) 

NUMBER 

OF 

LAYERS 

STANDARD 

DEVIATION 

(m/sec) 

SPREAD 

(km) 

0-3 

4 

11-15 

13-25 

3-19 

5 

5-8 

4-11 

19-34 

3 

10-15 

11-27 

34-99 

7 

15-20 

20-40 

99-199 

1 

3 

The  structure  is  best  resolved  in  the  3-19  km  range.  Both  the 
standard  deviations  and  the  spread  gradually  increase  below  this 
depth.  The  single  bottom  layer  is  well-resolved  because  of  the 
inclusion  of  long-period  phase  velocity  data,  but  the  spread  is  not 
defined  for  this  layer.  The  near  surface  layers  have  spread  greater 
than  depth,  an  indication  of  insufficient  information  content  in  the 
data  (in  this  case,  lack  of  high  frequency  data). 

The  frequency  range  that  was  obtained  from  the  synthetic 
seismograms  is  optimistic  for  actual  data.  A  typical  range  is 
0.02  Hz  -  0.12  Hz  (50-8  sec).  The  loss  of  long-period  data  reduces 
the  maximum  resolvable  depth,  and  such  a  restriction  would  probably 
limit  an  accurate  model  to  the  crust  and  very  top  of  the  mantle. 
However,  the  good  resolution  obtained  at  mid-crustal  depths  would  be 
retained. 

Another  problem  can  arise  because  of  the  +2nir  phase  ambiguity 
discussed  earlier.  The  result  of  choosing  an  incorrect  value  of  n 
will  cause  the  velocities  of  the  lower  layers  to  be  incorrect.  The 
velocities  of  the  upper  layers  are  more  strongly  constrained  by  the 
group  velocity  than  the  phase  velocity,  and  so  are  little  affected. 
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VI.  ATTENUATION,  MOMENT  ESTIMATES  AND  INVERSION 
FOR  Q  STRUCTURE 


Once  the  inverse  structure  has  been  determined,  synthetic 
spectral  amplitudes  can  be  computed  and  used  to  estimate  the  long 
period  source  level  u^oo  and  the  attenuation  coefficients  y(«).  The 
synthetic  test  seismograms,  of  course,  simulate  the  ideal  data  set 
and  therefore  produce  the  best  estimates  of  attenuation  and  source 
level  that  we  can  hope  to  obtain. 

For  the  synthetic  test,  we  used  inversion  model  2  (Figure  9) 
for  computation  of  inverse  synthetic  seismograms.  At  the  longest 
usable  period  of  200  seconds,  we  set  y(<*>)  *  0.  With  this 
constraint,  both  the  source  level  and  y(u>)  are  determined.  The 
estimated  value  of  differed  from  the  true  value  by  less  than  2 

percent.  A  comparison  of  the  "true"  y(w)  with  the  estimated  y(«»)  is 
shown  in  Table  1. 

In  order  to  assess  the  accuracy  of  the  attenuation 
coefficients,  we  inverted  them  (see  Equation  (4))  in  terms  of  a 
QB  (z)  model  (Figure  12).  The  inverted  Q  model  is  a  good 
approximation  to  the  "true"  Q  structure.  The  sharp  discontinuity  in 
Q  at  shallow  depths  has,  of  course,  been  smoothed  by  the  inversion 
procedure. 
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TABLE  1 


ACTUAL  AND  ESTIMATED  ATTENUATION  COEFFICIENTS 


Period  y  Estimated  (10~^/km)  y  True  (lO^/km) 


200 

0.00 

0.02 

100 

0.03 

0.05 

50 

0.15 

0.12 

40 

0.22 

0.18 

30 

0.35 

0.31 

25 

0.49 

0.47 

20 

0.77 

0.79 

15 

1.26 

1.41 

10 

2.47 

2.91 

5 

10.60 

11.70 
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jure  12. 


Inverted  Q-model  (bottom)  and  actual  Q-model  (top). 
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VII.  TEST  USING  HOGGAR  DATA 


Surface  wave  recordings  of  the  French  Sahara  event  SAPHIRE 
(27  February  1965)  were  used  to  test  the  S^  surface  wave  codes. 

The  recordings  were  from  two  WWSSN  stations,  Shirez,  Iran  (SHI)  and 

Addis  Ababa,  Ethiopia  (AAE).  These  seismograms  were  manually 

digitized  from  film-chip  reproductions  (Figure  13). 

These  two  seismograms  are  a  subset  of  a  larger  data  set 
processed  by  Rodi,  et  al.  (1978).  In  that  study,  dispersion  curves 
and  velocity  structures  were  derived.  However,  dispersion  curves 

were  found  by  narrow-band  filtering  rather  than  phase-matched 
filtering  and  no  attempt  was  made  to  perform  an  inversion  for 
Q  (z).  The  results  shown  in  this  report  were  found  independently 

D 

of  the  work  of  Rodi,  et  al.  (1978). 

The  SHI  seismogram  (Figure  13a)  shows  fairly  clean  dispersion, 
little  evidence  of  higher  modes,  and  moderate  background  noise.  A 
multiple  arrival  is  evident  at  about  1740  seconds.  Since  this 
arrival  is  not  clearly  separated  from  the  primary,  it  caused  a 
slight  perturbation  from  a  smooth  group  velocity  curve. 

The  AAE  seismogram  (Figure  13b)  shows  two  unusual  features: 
significant  inverse  dispersion  (high  frequencies  arrive  first),  and 
a  "boxcar"-l ike  amplitude  envelope.  Background  noise  and  possibly 
higher  modes  and  multi  paths  are  present. 

The  velocity  structures  inverted  from  the  two  seismograms  both 
show  a  smooth  gradient  from  surface  to  mantle  and  have  no  definite 
manifestations  of  a  crust/mantle  discontinuity  (Figure  14). 
However,  the  average  near-surface  velocity  for  the  SHI  path  is  much 
lower  than  that  for  the  AAE  path. 

The  comparison  of  observed  and  synthetic  spectral  amplitudes 
for  attenuation  estimates  led  to  several  difficulties.  Most 
notably,  the  AAE  spectrum  peaked  at  intermediate  periods,  and  so 
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constraining  y(«»)  *  0  at  the  longest  period  would  have  resulted  in 
negative  apparent  attenuation  at  shorter  periods.  We  then  changed 
our  constraint  so  that  y(«)  is  positive  for  all  frequencies. 

The  resulting  Q  (z )  structures  are  shown  in  Figure  15.  The 

D 

models  are  quite  different  in  spite  of  the  fact  that  the  paths 
differ  in  azimuth  by  only  35  degrees.  Needless  to  say,  we  can  only 
express  low  confidence  in  these  results.  The  observed/synthetic 
spectral  amplitude  ratio  showed  differences  that  are  too  great  to  be 
attributed  to  simple  attenuation,  and  the  inverted  Q  ( z ) 

O 

structures  are  implausible.  Nevertheless,  the  synthetic  seismograms 

resulting  from  the  inverted  0 ( z)  and  Q.(z)  structures  generally 

s 

match  the  actual  seismograms  quite  well  (Figure  16).  The  peak 
amplitude  for  SHI  is  about  one  cycle  too  late,  probably  due  to 

difficulties  in  measuring  the  dispersion  curve  near  the  Airy  phase. 
Note  tnat  the  suspected  multiple  arrival  mentioned  earlier  is  not 
present  in  the  synthetic. 

Background  noise  is  likely  to  be  a  major  contributing  factor 
in  the  observed  errors,  particularly  since  the  data  was 
hand-digitized.  However,  we  suspect  that  an  even  more  serious 
problem  arises  from  lateral  variations  in  structure  along  the  path. 
No  reasonable  plane-layered  model  would  accurately  replicate  the 

observed  spectrum.  It  is  quite  likely  that  a  simple  plane-layered 
model  is  inadequate  when  the  path  involves  complex  tectonic 

features.  In  our  examples,  station  SHI  is  near  the  Zagros  Suture 

Zone,  while  station  AAE  is  in  close  proximity  to  the  East  Africa 
Rift  Zone. 

Our  current  work  is  focused  on  identifying  such  problems  and 
implementing  appropriate  Improvements  to  the  basic  algorithms. 
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Inverted  Q0(z)  structures  from  SAPHIRE  seismograms: 
(b)  AAE.  8 
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VIII.  COMPUTER  CODES 


3 

The  S  surface  wave  processing  is  broken  down  into  three 
separate  programs:  TELVEL,  for  computing  dispersion  curves  and 

spectral  amplitudes;  SWIP,  for  inverting  dispersion  curves  and  y(m) 
for  b(z)  and  Q _ ( z) ;  and  SURWAV,  for  generating  synthetic 

D 

seismograms  (Figure  2). 

TELVEL  is  an  interactive  program  with  graphics  designed  to 

work  on  a  Tektronix  4014  CRT  storage  terminal.  Operator  interaction 

is  necessary  to  recognize  multipath  arrivals  and  to  select 

appropriate  time  windows.  Currently,  TELVEL  comprises  about  30 

subroutines  and  requires  65  K  core  storage.  Approximately  20  K  of 

the  storage  requirement  is  due  to  use  of  the  proprietary  graphics 
R  3 

package  DISSPLA  .  S  is  now  developing  its  own  site* independent 
graphics  package  to  enhance  program  transportability  and 

simultaneously  reduce  storage  requirements  (to  about  10  K  for 
graphics). 

SWIP  is  currently  a  batch-run  program  with  input  parameters 
specified  in  namelist  statements  and  free-format  reads.  It  is  not 
configured  to  directly  read  the  TELVEL  output  file,  but  instead 
requires  text  editing  by  the  operator.  SWIP  writes  results  to  a 
disk  file. 

SURWAV  is  also  currently  a  batch-run  program  with  namelist 
inputs.  It  writes  its  synthetic  seismograms  in  SDL  format. 

The  three  programs  mentioned  above  were  developed 
independently  of  each  other;  the  interfacing  and  data  transfer 
between  the  programs  is  currently  being  streamlined.  File  structure 
is  being  reformatted  to  be  consistent  between  the  programs  and  to 
minimize  operator  effort. 

The  two  batch  codes  SWIP  and  SURWAV  are  also  being  redesigned 
to  be  interactive  and  to  automate  as  many  steps  as  possible.  While 
this  requires  a  significant  amount  of  code  development,  it  should 
greatly  reduce  the  amount  of  time  necessary  to  obtain  path 
corrections  after  the  modifications  are  completed. 
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IX.  OBJECTIVES  AND  RECOMMENDATIONS 


Our  primary  objective  under  this  contract  is  to  develop  and 
transfer  to  VSC  the  software  required  to  obtain  path  corrections. 
We  are  now  in  the  process  of  modifying  our  existing  codes  to  make 
them  easy  to  use  and  to  make  the  determination  of  path  corrections 
as  rapid  as  possible.  As  soon  as  this  is  accomplished,  the  codes 
will  be  transferred  to  VSC  and  adapted  to  be  consistent  with  the  VSC 
file  structure  and  computer  facilities.  They  will  also  be  adapted 
to  interface  with  other  codes  at  VSC  such  as  the  moment  tensor 
inversion  package  being  supplied  by  Sierra  Geophysics  which  will  use 
the  path  corrections. 

We  are  also  attempting  to  improve  the  method  of  estimating  the 
attenuation  coefficients.  In  general,  the  Rayleigh  wave  spectra  do 
not  contain  enough  information  to  determine  the  attenuation 
coefficients  frequency  by  frequency.  It  is  necessary  to  restrict 
the  allowable  Q  models  in  some  way.  There  are  several  possible 
methods  of  doing  thisi  If  detailed  attenuation  studies  have  been 
done  in,  or  near,  the  area  of  interest,  a  published  Q  model  may  be 
available.  Another  possibility  is  to  make  Q  a  function  of  shear 
velocity  or  shear  modulus  since  low  Q  zones  typically  correlate  with 
regions  of  lower  shear  velocity.  After  an  initial  trial  Q  model  is 
selected,  it  may  be  modified  by  changing  a  restricted  number  of  free 
parameters  and  comparing  the  synthetic  with  the  observed  seismograms 
to  obtain  an  improved  Q  model  and  improved  attenuation  coefficients. 

In  addition  to  developing  and  transferring  the  codes,  we  are 
preparing  user  manuals  to  accompany  them.  While  the  codes  are 
designed  to  obtain  path  corrections,  they  are  useful  for  a  wide 
variety  of  other  tasks  and  should  be  a  valuable  addition  to  the  VSC 
computing  facility. 
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